'''
Created on Sep 18, 2010

@author: oabalbin
'''


'''
Create a dictionary of intervals: key; interval hg18_chr_start_end_strand = name
    using the bed file

read the sequence file save using the Bioseq fasta 
seq name, join the sequences using the name of sequence. 
    Remember to index the fasta file using the dict of Bioseq

'''

import numpy as np
from Bio.Seq import Seq
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from optparse import OptionParser

from collections import deque, defaultdict

def create_list_of_intervals(intervals_bed_file):
    
    intervals_bed_file = open(intervals_bed_file)
    interval_list=defaultdict(deque)
    for line in intervals_bed_file:
        fields = line.strip('\n').split('\t')
        name, score, strand=fields[3],fields[4],default_strand(fields[5])
        dname='hg18_'+fields[0]+'_'+fields[1]+'_'+fields[2]+'_'+strand
        interval_list[name].append(dname)
        
    intervals_bed_file.close()
    return interval_list

def default_strand(x):
    if not x:
        return '+'
    elif x=='.':
        return '+'
    else:
        return x

def combine_exons_sequence(inputfile_name, gene_list, outfile_name):
    
    outfile = open(outfile_name,'w')
    
    record_dict = SeqIO.index(inputfile_name, "fasta")
    #c=0
    for gname, exons_list in gene_list.iteritems():
        gene_sequence=''
        geneloc2=''
        geneloc=[]
        for exon in exons_list:
            geneloc2 += ",".join(exon.split('_')[1:4]).replace(',','_')+'~' 
            geneloc.append(exon.split('_')[1:4])
            gene_sequence += record_dict[exon].seq
        
        outfile.write('>'+gname+'_'+",".join(map(str,geneloc[0][0:2])).replace(',','_')+\
                      '_'+str(geneloc[-1][2])+\
                      '_'+','.join(map(str,exon.split('_')[-1]))+'\n')

        i,ft,seqlen=0,50,len(gene_sequence)
        seqtail=np.mod(seqlen,ft)
            
        while (i < (seqlen-seqtail)):
            outfile.write(str(gene_sequence[i:i+ft])+'\n')
            i+=ft
        if seqtail>0:
            outfile.write(str(gene_sequence[i:i+seqtail])+'\n')
        
        '''                
        c+=1
        if c>10:
            break
        '''    
    outfile.close()
    

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-t", "--transcripts", dest="transcripts",
                            help="Bed file with the location of transcripts of interest")
    optionparser.add_option("-s", "--sequences", dest="sequences",
                            help="fasta file with sequences for each transcript")
    optionparser.add_option("-o", "--transcriptome", dest="transcriptome",
                            help="fasta file with the join sequences of the trasncripts")

    (options, args) = optionparser.parse_args()

    '''
    bed_file = '/exds/users/oabalbin/projects/mirnas/test/test1/2010_09_01_14_30_classifier_transcripts_exons.bed'
    inputfile_name = '/exds/users/oabalbin/projects/mirnas/test/test1/alu_exons_genenomic_sequences.fa'
    outfile_name = '/exds/users/oabalbin/projects/mirnas/test/test1/alu_exons_transcriptome_sequences.fa'
    '''
    bed_file = options.transcripts
    inputfile_name = options.sequences
    outfile_name = options.transcriptome
    gene_list = create_list_of_intervals(bed_file)
    combine_exons_sequence(inputfile_name, gene_list, outfile_name)
